function [A,S,dx,dy] = twoDcoef( dat, xcm, xfm, ycm, yfm, mt, src, BC )
% function A = twoDcoef( dat, xcm, xfm, ycm, yfm, mt, src, BC )
%   This function produces the coefficient matrix and source for solution of the
%   two-dimensional diffusion equation.  Currently, one-group.

% expand the bc's
BCL = BC(1);
BCR = BC(2);
BCB = BC(3);
BCT = BC(4);

% define data components
dc = dat(:,1); % diffusion coefficient
sa = dat(:,2); % absorption cross-section

% number of fine meshes in x and y coordinates
N = sum(xfm); 
M = sum(yfm);
% number of coarse meshes in x and y coordinates
CN = length(xfm);
CM = length(yfm);
% compute dx and dy vectors
dx = zeros(N,1);
dy = zeros(M,1);
j = 0;
for i = 1:CN
    dx( (j+1):(j+xfm(i)) ) = (xcm(i+1)-xcm(i))/xfm(i);
    j = sum(xfm(1:i));
end
j = 0;
for i = 1:CM
    dy( (j+1):(j+yfm(i)) ) = (ycm(i+1)-ycm(i))/yfm(i);
    j = sum(yfm(1:i));
end

lenk = (N+1)*(M+1);
AC = zeros( lenk, 1 );
AL = zeros( lenk-1, 1 );
AR = AL;
AB = zeros( lenk-N-1, 1 );
AT = AB;
S = zeros( (N+1)*(M+1), 1 );

% interior coefficients
for i = 2:N 
    for j = 2:M
        for ci = 1:CN
            if i-1 <= sum(xfm(1:ci)), mtx = ci; break; end
        end
        for ci = 1:CM
            if j-1 <= sum(yfm(1:ci)), mty = ci; break; end
        end
        for ci = 1:CN
            if i  <= sum(xfm(1:ci)), mtxo = ci; break; end
        end
        for ci = 1:CM
            if j  <= sum(yfm(1:ci)), mtyu = ci; break; end
        end        
        mtt  = mt(mtx,mty);   mtu  = mt(mtx,mtyu);
        mto  = mt(mtxo,mty);  mtb  = mt(mtxo,mtyu);
        
        k = i+(j-1)*( N+1 ); % 10 being the number of x ints
        alft = -0.5*( dc(mtt)*dy(j-1) + dc(mtu)*dy(j) ) / dx(i-1);
        argt = -0.5*( dc(mto)*dy(j-1) + dc(mtb)*dy(j) ) / dx(i);
        abot = -0.5*( dc(mtt)*dx(i-1) + dc(mto)*dx(i) ) / dy(j-1);
        atop = -0.5*( dc(mtu)*dx(i-1) + dc(mtb)*dx(i) ) / dy(j);    
        a    = 0.25*( dx(i-1)*dy(j-1)*sa(mtt) + dx(i)*dy(j-1)*sa(mto) + ...
               dx(i-1)*dy(j)*sa(mtu) + dx(i)*dy(j-1)*sa(mtb) ) - (alft+argt+abot+atop);
        AL( k-1 )   = alft;  
        AR( k )     = argt;   
        AB( k-N-1 ) = abot;    
        AT( k )     = atop;
        AC( k )     = a;
        S(k)        = 0.25*( dx(i-1)*dy(j-1)*src(mtx,mty) + dx(i)*dy(j-1)*src(mtxo,mty) + ...
                             dx(i-1)*dy(j)*src(mtxo,mty) + dx(i)*dy(j)*src(mtxo,mtyu) );              
    end
end

% left and right edges
for j = 2:M
    for ci = 1:CN
        if 2 <= sum(xfm(1:ci)), mtx = ci; break; end
    end
    for ci = 1:CM
        if j-1 <= sum(yfm(1:ci)), mty = ci; break; end
    end
    for ci = 1:CN
        if N-1  <= sum(xfm(1:ci)), mtxo = ci; break; end
    end
    for ci = 1:CM
        if j  <= sum(yfm(1:ci)), mtyu = ci; break; end
    end       
    i = 1;
    k = i+(j-1)*( N+1 );
    if BCL == 1
        AR( k ) = -1;             
        AC( k ) = 1;  
    else
        % average d in vert
        AR( k ) = -0.25 * ( dc(mt(1,mty)) + dc(mt(1,mtyu)) ) / dx(i);    
        AC( k ) = 0.25 + 0.5 / dx(i) * ( dx(i)*dy(j-1)*dc(mt(1,mty)) + dx(i+1)*dy(j-1)*dc(mt(mtx,mty)) +...
                                       dx(i)*dy(j)*dc(mt(1,mtyu))  + dx(i+1)*dy(j)*dc(mt(1,mtyu)) ) / ...
                                        ( dx(i)*dy(j-1)+dx(i+1)*dy(j-1)+dx(i)*dy(j)+dx(i+1)*dy(j) );  
    end
    i = N+1;
    k = i+(j-1)*( N+1 );
    if BCR == 1
        AL( k-1 ) = -1;              
        AC( k )   = 1;
    else
        AL( k-1 ) = -0.25 * ( dc(mt(end,mty)) + dc(mt(end,mtyu)) ) / dx(N);
        AC( k ) = 0.25 + 0.5 / dx(N) * ( dx(N)*dy(j-1)*dc(mt(end,mty)) + dx(N-1)*dy(j-1)*dc(mt(mtxo,mty))  +...
                                       dx(N)*dy(j)*dc(mt(end,mtyu))  + dx(N-1)*dy(j)*dc(mt(mtxo,mtyu)) ) / ...
                                        ( dx(N-1)*dy(j-1)+dx(N)*dy(j-1)+dx(N-1)*dy(j)+dx(N)*dy(j) );      
    end
end

% top and bottom edges
for i = 2:N
    for ci = 1:CN
        if i-1 <= sum(xfm(1:ci)), mtx = ci; break; end
    end
    for ci = 1:CM
        if 2 <= sum(yfm(1:ci)), mty = ci; break; end
    end
    for ci = 1:CN
        if i  <= sum(xfm(1:ci)), mtxo = ci; break; end
    end
    for ci = 1:CM
        if M-1  <= sum(yfm(1:ci)), mtyu = ci; break; end
    end 
    j = 1;
    k = i+(j-1)*( N+1 );  
    if BCB == 1
        AT( k )     = -1;             
        AC( k  )    = 1;  
    else
        AT( k )     = -0.25 * ( dc(mt(mtx,1)) + dc(mt(mtxo,1)) ) / dy(1);   
        AC( k )     = 0.25 + 0.5 / dy(1) * ( dx(i-1)*dy(j)*dc(mt(mtx,1)) + dx(i)*dy(j)*dc(mt(mtxo,1)) +...
                                       dx(i-1)*dy(j+1)*dc(mt(mtx,mty))  + dx(i)*dy(j+1)*dc(mt(mtxo,mty)) ) / ...
                                        ( dx(i-1)*dy(j)+dx(i)*dy(j)+dx(i-1)*dy(j)+dx(i)*dy(j+1) );            
    end 
    j = M+1;
    k = i+(j-1)*( N+1 );
    if BCT == 1
        AB( k-N-1 ) = -1;            
        AC( k )     = 1; 
    else
        AB( k-N-1 ) = -0.25 * ( dc(mt(mtx,end)) + dc(mt(mtxo,end)) ) / dy(M);    
        AC( k )     = 0.25 + 0.5 / dy(M) * ( dx(i-1)*dy(M-1)*dc(mt(mtx,mtyu)) + dx(i)*dy(M-1)*dc(mt(mtxo,mtyu)) +...
                                       dx(i-1)*dy(M)*dc(mt(mtx,end))  + dx(i)*dy(M)*dc(mt(mtxo,end)) ) / ...
                                        ( dx(i-1)*dy(M-1)+dx(i)*dy(M)+dx(i-1)*dy(M-1)+dx(i)*dy(M) );         
    end
end

% finally, the corners
for ci = 1:CN
    if i-1 <= sum(xfm(1:ci)), mtxL = ci; break; end  % picks out material for xfm(2)
end
for ci = 1:CM
    if 2 <= sum(yfm(1:ci)), mtyB = ci; break; end % picks out material for yfm(2)
end
for ci = 1:CN
    if i  <= sum(xfm(1:ci)), mtxR = ci; break; end % picks out material for xfm(end-1)
end
for ci = 1:CM
    if M-1  <= sum(yfm(1:ci)), mtyT = ci; break; end % picks out material for yfm(2)
end 
i=1; j=1; k = i+(j-1)*( N+1 );      % BOTTOM LEFT
if (BCL==1&&BCB==1)
    AR( k ) = -1; %right  
    AT( k ) = -1; %top
    AC( k ) = 2;  %center
elseif(BCL==1&&BCB==0)
    AR( k ) = -1; %right  
    AC( k ) = 1; %center   
elseif (BCL==0&&BCB==1)
    AT( k ) = -1; %top
    AC( k ) = 1; %center   
else    
    AR( k ) = -0.25 / dx(1) * ( dc(mt(1,1)) + dc(mt(mtxL,1)) ) ; %right  
    AT( k ) = -0.25 / dy(1) * ( dc(mt(1,1)) + dc(mt(1,mtyB)) ); %top
    AC( k ) = 0.5 + 0.5*( dc(mt(1,1))/dx(1)+dc(mt(1,1))/dy(1) ); %center  
end
i=1; j=M+1; k = i+(j-1)*( N+1 );    % TOP LEFT
if (BCL==1&&BCT==1)
    AR( k )     = -1; %right  
    AB( k-N-1 ) = -1; %bot
    AC( k )     = 2;  %center
elseif(BCL==1&&BCT==0)
    AR( k ) = -1; %right  
    AC( k ) = 1;  %center   
elseif (BCL==0&&BCT==1)
    AB( k-N-1 ) = -1; %bot
    AC( k )     = 1;  %center   
else    
    AR( k )     = -0.25 / dx(1) * ( dc(mt(1,end)) + dc(mt(mtxL,end)) ); %right  
    AB( k-N-1 ) = -0.25 / dy(end) * ( dc(mt(1,end)) + dc(mt(1,mtyT)) ); %bot
    AC( k )     = 0.5 + 0.5*( dc(mt(1,end))/dx(1) + dc(mt(1,end))/dy(end) ); %center  
end
i=N+1; j=1;  k = i+(j-1)*( N+1 );   % BOTTOM RIGHT
if (BCR==1&&BCB==1)
    AL( k-1 ) = -1; %left  
    AT( k )   = -1; %top
    AC( k )   = 2;  %center
elseif(BCR==1&&BCB==0)
    AL( k-1 ) = -1; %left  
    AC( k )   = 1; %center   
elseif (BCR==0&&BCB==1)
    AT( k )     = -1; %top
    AC( k )                   = 1; %center   
else    
    AL( k-1 ) = -0.25 / dx(end) * ( dc(mt(end,1)) + dc(mt(mtxR,1)) ) ; %left  
    AT( k )   = -0.25 / dy(1) * ( dc(mt(end,1)) + dc(mt(end,mtyB)) ); %top
    AC( k )   = 0.5 + 0.5*( dc(mt(end,1))/dx(end)+dc(mt(end,1))/dy(1) ); %center  
end
i=N+1; j=M+1;  k = i+(j-1)*( N+1 );   % TOP RIGHT
if (BCR==1&&BCT==1)
    AL( k-1 )   = -1; %right  
    AB( k-N-1 ) = -1; %bot
    AC( k )     = 2;  %center
elseif(BCL==1&&BCT==0)
    AL( k-1 ) = -1; %right  
    AC( k )   = 1; %center   
elseif (BCL==0&&BCT==1)
    AB( k-N-1 ) = -1; %bot
    AC( k )     = 1; %center   
else    
    AL( k-1 )   = -0.25 / dx(end) * ( dc(mt(end,end)) + dc(mt(mtxL,end)) ); %right  
    AB( k-N-1 ) = -0.25 / dy(end) * ( dc(mt(end,end)) + dc(mt(end,mtyT)) ); %bot
    AC( k )     = 0.5 + 0.5*( dc(mt(end,end))/dx(1) + dc(mt(end,end))/dy(end) ); %center  
end    

%---establish sparse matrix for solution
% pad end of AB,AL
AB = [AB' zeros(1,N+1)]';
AL = [AL' 0]';
% pad top of AR,AT
AT = [zeros(1,N+1) AT']';
AR = [0 AR']';
% form sparse matrix for solution
A = spdiags([AB AL AC AR AT],[-N-1 -1 0 1 N+1],(N+1)*(M+1),(N+1)*(M+1));